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Abstract A computational method for magnetically guided drug delivery is pre¬ 
sented and the results are compared for the aggregation process of magnetic par¬ 
ticles within a fluid environment. The model is developed for the simulation of 
the aggregation patterns of magnetic nanoparticles under the influence of MRI 
magnetic coils. A novel approach for the calculation of the drag coefficient of ag¬ 
gregates is presented. The comparison against experimental and numerical results 
from the literature is showed that the proposed method predicts well the aggrega¬ 
tions in respect to their size and pattern dependance, on the concentration and the 
strength of the magnetic held, as well as their velocity when particles are driven 
through the fluid by magnetic gradients. 

Keywords MRI guided drug delivery, Aggregations, Magnetic nanocapsules 


1 INTRODUCTION 

Erom the beginnings of 1970s researchers were studying the concept of magnetic 
guided drug delivery method [^[^. The concept of this method is to attach the 
drug to the micro- or nanoparticles and then to inject them to the bloodstream. Eor 
the guidance to the targeted area, a Magnetic Resonance Imaging (MRI) device is 
needed. By making use of the magnetic guided drug delivery method, the quantity 
of the drug required to reach therapeutic levels is being reduced. Also, the drug 
concentration at targeted sites is increased. 

The efficiency of this concept depends on the materials of the particles (cobalt 
and manganese ferrites, encapsulated or coated [^), the blood flow rate and the 
intensity of the magnetic held [^. It is found that, this method is more efficient 
in small blood vessels with low blood flow rates than large blood vessels, where 
the flow rate is higher. A small size of particles implies a magnetic response of 
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reduced strength and as a result, it is difficult to drive particles and keep them 
in the targeted region [^. In order to overcome this difficulty the use of magnetic 
particles to create aggregations is proposed in order to increase the magnetic re¬ 
sponse [^. The total magnetic moment of clusters is higher than that of isolated 
particles, and therefore, cluster are more magnetically responsive [^. The ratio of 
the velocity of the aggregates to the velocity of an individual particle was found to 
reach a constant value independent of the aspect ratio value . The magnetization 
of nanospheres increases as a result of the increase in the concentration of matter 
along the direction of the field. This increase leads to the observed change in veloc- 
ity [^. When the aggregations reach the targeted region, they break up into single 
particles. This break up can be achieved by using superparamagnetic particles that 
lose the magnetization after moving out from the magnetic field [^. The size of 
aggregates is a very important parameter, since large aggregations could form clots 
in small arteries. On the other hand, a small aggregation would be dragged away 
from the cardiovascular system circulation. Therefore, the size of aggregation has 
an important role in the effort to create an efficient propulsion through the blood 
vessels. Magnetic drug targeting is possible in treating superficially located dis¬ 
eases, such as tumors and vascular diseases. The treatment of tissues located away 
from the body surface is not possible because the magnetic force decrease with the 
increasing distance from the electromagnetic coil [^. In order to overcome this 
difficulty magnets are implemented near the targeted sites 0 0 


Since an analytical study of aggregation process is impossible to be developed, 
a numerical model for magnetically guided drug delivery is attended. The method 
that is proposed here simulates the movement of aggregated magnetic particles in 
a fiuid environment. The numerical model can simulate the number of resulting 
aggregations whose size and pattern depend on the concentration and the strength 
of the magnetic field. The forces acting on a particle are described in detail in 
Section 2. The model is validated through a comparison with two benchmark test 
cases of Ref. and Ref. and the results are discussed in Section 3. Finally, 
discussion and conclusions are presented in Sections 4 and 5, respectively. 


2 METHODS 

For the propulsion model of the particles, six major forces are considered, i.e. the 
magnetic force from MRIs Main Magnet static field as well as the Magnetic field 
gradient force from the special Propulsion Gradient Coils. The static field caters 
for the aggregation of nanoparticles while the magnetic gradient navigates the 
agglomerations. Moreover, the contact forces among the aggregated nanoparticles 
and the wall, and the Stokes drag force for each particle are considered, while only 
spherical particles are used here. Finally, gravitational forces due to gravity and 
the force due to buoyancy are added. 

The motion of particles are given by the Newton equations: 

du ■ 

~ F'magA T ^ncA T F'tcA T ^ drag A T ^boyA T i (1) 


li 


dujj 

dt 


= M 


drag A 


+ ^ con A + T 


mag A 


( 2 ) 
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where, the index i stands for the particle i. The bold variables are vectors. The 
velocity is Ui and the rotational velocity is The mass of particle z is t stands 
for the time and the mass moment of inertia matrix is The linear and angular 
acceleration are and respectively. The total magnetic force is Fmag.i- 
Fnc_i and Ffcj are the normal and tangential contact forces, respectively. F^^agA 
stands for the hydrodynamic drag force. Fi^^yj and Wi are the buoyancy and the 
weight forces. M^^^agA and M^onA stand for the drag and the contact moments, 
respectively. T^agA is the torque due to the magnetic field at the position of 
particle i. 

In the following section, the major forces that are taking into account during 
the simulation are described in detail. 


2.1 Forces acting on particles 


1. Magnetic Forces: Magnetic forces, F^agA^ exerted on the i particle are given 

by [igi^ : 


magji 


= F 


intmagA 


+ F. 


ismagjL 


(3) 


where, FintmagA is the magnetic force due to the interaction of the zth particle 
with the magnetic field, and FismagA is the magnetic force acting on the zth 
particle due to its interaction with the surrounding magnetic particles. 

The force FintmagA in liie volume of the ferromagnetic body, V, is given by: 


FintmagA — ^ 

The force FismagA is given by: 

N 

FismagA — ^ ^ Figmag-ji 
3 

The numerical model for the magnetic forces is given in [13| . 

2. Fluid Forces: Each particle is subjected to the drag force that reads as: 

1 2 

^ drag A ~ 2^^ 


(4) 

( 5 ) 


(6) 


where, p is the density of the fluid, u is the speed of the particle i relative to 
the fluid and A is the reference area which equals to zrr^, where r stands for 
the radius of spherical particle. is the drag coefficient given by: 




24 [1 + 0.15i?e°-®^^] 
Re 


(7) 


where. Re is the Reynolds number based on the particle diameter [17| . 

In the present work, a new model is developed for the substitution of the 
reference area A, in Eq. (6), with the effective area A^ff. This substitution takes 
into account the fact that the downstream spheres of the chain lie partially in 
the wake of the upstream one, therefore the reference area A in Eq. (6) of each 
sphere must be reduced to each real area that is exposed to the flow. On this 
basis the area A^ff can be calculated as follows : 
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Fig. 1: Relevant velocity and drag force of a particle. 


(a) Sweep all spheres. 

(b) Find the tangent sphere to the current one. 

(c) A straight chain assumption considered where the sphere (1) is tangent 

only to two spheres namely the upwind and downwind. Sphere (2) which 
lies upstream to sphere (1) is found out by checking the inner product 
32 * (Uobj ~ ^ fluid) fo be positive, see Figs. Ill and Where, Uobj stands 

for the velocity of the sphere and Ustands for the velocity of the fluid 

(d) Project the area of spheres (2 and 1) on the plane \^obj ~ ^ fluid that is 
cross sphere (2) at each center, see Fig. It- This yields into two intersected 
circles, as shown in Fig. 

(e) Change to plane coordinate system (2D). 

(f) Calculate the overlapping area (white domain in Fig.[^. The mathematical 
formula reads as: 

^overlapping — domain(2?>A) + do?nam(134) — F^(1234) (8) 

(g) Substitute Eq. (6) with the equation: 

^drag-i 2^^ ^d^eff (9) 

In this way, the real reference area A^ff of sphere exposed to the flow, i.e.: 

Aeff ~ ~ Eoverlapping (10) 

is taken into account. 

3. Collision Forces: Each particle interacts with other particles through contact. 
In our model we use the Discrete Element Method (DEM) to calculate the 
motion of particles [^. DEM is a numerical model capable of describing the 
mechanical behaviour of assemblies of spheres and computing the motion and 
shear effect of a large number of small particles. In DEM, particles are ap¬ 
proximated as rigid bodies and the interactions between them are explicitly 
considered [^. 

4. Body Forces: Both gravity and buoyancy are being included in the calculation 
of the body force. The force is addressed by 

4 3 

F'grav-i T ^boy-i 3^^* Pf)^ (H) 
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(a) (b) 

Fig. 2: a) 3 (Xs,Ys, Zs) stands for the projection of sphere center 2 (X 2 , 12 ,^ 2 )) 
on the plane Uobj ~ ^fluid crossing sphere 1 at each center (Xi,li,Zi), and b) 
1 ' {Xif ,Yif, Zif) stands for the projection of sphere center 1 (Xi,Yi, Zi)) on the 
plane Uobj — ^fluid crossing sphere 2 at each center (X 2 , 12 ,^ 2 )- 



Fig. 3: Circular sector of two spheres Eoverlapping = domam(234) domain (ISA) — 
F;(1234). 


where, pi and pf are the density of the particle i and the fluid, respectively, 
is the radius of the particle i and g is the acceleration due to gravity. 


2.2 Magnetic field and interaction domain 

The magnetic field B in the MRI bore is given by 

B = Bo + G + Bi (12) 

where, Bq is the MRI superconducting magnet field that is constant and uniform. 


20 


G is the gradient field and Si is the time dependent radio frequency field 

It is known that, the magnetic interaction force, F, in the parallel direction 
between two spheres is inversely proportional to the fourth power of the separation 
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Relevant distance of particles (m) 


Fig. 4: Velocity of particles when approach under the force of a 1.5 T constant 
horizontal magnetic field from initial distance of 15 particle diameters of 11 jim. 


distance 21l22 


Mi Mi 

F oc 


(/l + 0-2 + 


(13) 


where, Mi and Mj are the magnetic moments of the center of each sphere and h 
is the nearest distance between the surfaces of two spheres with radii and aj . 

Although, a weak interaction is being observed for particles that are more 
than five radii apart, as time goes by, the interaction force is getting stronger, 
because the particles are getting faster close to each other, as shown in Fig. 
for the velocity between two approaching particles. Each particle interacts with 
all the other which are in the same area. The magnetic interaction domain is 
defined by the area of the domain that the particles are located in. This is a time 
consuming method, because every time the magnetic moments of all particles are 
calculated. In order to accelerate the method, the magnetic moments that are 
exerted on particle i are calculated from particles that are located in distance of 
15 diameters. Then, only the magnetic moments that are located in distance of 
10 diameters are stored. In this way, the magnetic moments of particles that are 
located in distance of 10 diameters are automatically recalled in order to be used 
for the calculation of the magnetic moments of other particles. Simulations of 64 
and 125 particles of the method which was described above show agreement in half 
time with the method of calculating the magnetic moments exerted on particle i 
from all the other particles. 


2.3 Numerical Method 

The OpenFoam platform was used for the calculation of the fiow field and the 
uncoupled equations of particles motion |^. The simulation process reads as 
follow: Initially, the fiuid fiow is found using the incompressible Navier-Stokes 
equations and the pressure correction method. Upon finding the fiow field (pressure 
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Table 1: Simulation parameters 


Flow domain and computational grid of Case 1 (3D) 

No 

Concentr. 

Particles 

Volume 

Dimensions in (m) X 10^ 

Number of 


{mg/ml) 


(m3) 

oi X X y X z directions 

cells (x,y,z) 

1 

0.563 

311 

4.185 X 10“ 

7.48 

7.48 

7.48 

34 X 34 X 34 

2 

1.125 

385 

2.596 X 

6.38 

6.38 

6.38 

29 X 29 X 29 

3 

2.25 

336 

1.133 X 

4.84 

4.84 

4.84 

22 X 22 X 22 

4 

4.5 

307 

5.174 X 10“^^ 

3.96 

3.96 

3.30 

18 X 18 X 15 

Flow domain and computational grid of Case 2 (2D) 

5 1 

25 

1 55 

1 2.038 X 10-^^ 

1 1.924 

1 1.924 1 

0.055 1 

7x7 


and velocity) the motion of particles is evaluated by the Lagrangian method by 
solving Eqs. (1) and (2) along the trajectory of each particle. The equations are 
solved in time by the Euler time marching method. The stability of the algorithm 
is guaranteed through a time step of the order of 

The present numerical methods were validated against the results from Refs. 

In order to perform the comparisons, two series of simulations with the fol¬ 
lowing computational domains and grid distributions were selected: a) In the first 
case, four water solutions with different concentrations (0.563, 1.125, 2.25 and 4.5 
in mg/ml) and by using about 300 — 400 particles were simulated under the mag¬ 
netization of a uniform magnetic field of 0.4 T in a stationary fluid. The spacing 
of the 3D computational grid of this case was 2 x diameter of the particle in each 
direction, b) In the second case, a solution with concentration of 25 mg/ml of 
55 polystyrene magnetic particles with diameter of 5.5 gm was simulated in dis¬ 
tilled stationary water under a constant magnetic field of 0.005T and a gradient 
magnetic field of (5 = 1.4 T/m along the a:— axis. The spacing of the 2D domain 
was 5 X diameter of the particle in each direction. The summary of the domain 
parameters for the first and the second case is shown in Table 


3 RESULTS 

Two series of simulations were performed according to the parameters that are 
described in the previous section. In the first case, which is initially addressed 
by Mathieu and Martel [^, the aggregation of Fes04 particles dispersed in a 
stationary distilled water is investigated under the influence of a constant magnetic 
field. The density of the particle and the fluidic environment was 1087 kg/w? and 
1000 kg/vr? ^ respectively, while the diameter of the particles was considered to be 
11 gm. The relative magnetic permeability of the Fe^O^ particles was 12.3 and 
the magnetic permeability of the fluid was An x 10 “^ A/m. The Young’s modulus 
of the material was 10 ^ N/m?^ the tangential stiffness was 10 Nsm~^ and the 
coefficient of friction and the Poisson ratio was 0.5 for the two properties. 

Eour water solutions with particle concentrations of 0.563 mg/ml^ 1.125 mg/ml^ 
2.25 mg/ml and 4.5 mg/ml were simulated under a uniform transverse magnetic 
field of magnitude Bq = 0.4 T. No external magnetic gradient is applied during this 
simulation and the initial positions of the microparticles are randomly generated. 
At the beginning of each simulation, for t = Os, the magnetic field is Rq = 0 T 
and under the influence of the uniform magnetic field the particles (magnetic 
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0.009 
0.008 
0.007 
0.006 
0.005 
^ 0.004 
0.003 
0.002 
0.001 


Fig. 5: Green (lighter) particles are shown the initial particle position at t = Os, 
while red (darker) particles correspond to t = 5 ms where the particles chains are 
fully formed parallel to the steady magnetic held. 



dipoles) are formed into chains that are oriented parallel to the magnetic held 
within t = 5 ms, as it is shown in Fig. 

Due to the external magnetic held, the magnetic particles are acting like mag¬ 
nets and the particle motion is attributed to gravity and magnetic interactions 
forces, since no huid how exists. The size of the aggregations that are formed as 
a function of particle numbers and concentrations is presented in Table Results 
from concentration 1.125 mgjml, for example, show that four aggregates of length 
10 particles may be formed due to the application of this magnetic held. It is 
observed that the size of the aggregations is increased as concentration increases 
and for 4.5 mg/ml aggregates of 35 particles may be found. The results from the 
simulation are compared against the experimental data of Ref. as it is depicted 
in Fig. and tabulated in Table for the mean length of the aggregates and their 
standard size deviation. 

Figure shows the length of aggregates major axis whose average value fol¬ 
lows an upward trend as a function of suspension concentration. The major axis 
length corresponds to the longest dimension of the aggregates measured. Larger 
aggregates are detected in increasing numbers as the concentration of the sus¬ 
pension increases. For the concentration of 0.563 mg/ml the numerical results are 
very close to the experimental, while for the concentration of 1.125 mg/ml results 
are similar compared to the experimental data in terms of mean length. For the 
higher concentrations an underestimation may be observed between the results of 
simulations and the experimental one. This difference is due to the small numbers 
of particles (up to 380 particles) that were used in the simulation in comparison 
in contrast to the measurements that were conducted using a number of particles 
in the order of 10^. 

The second case that is studied here, which is initially addressed by Vartholomeos 
and Mavroidis [^, investigates the flow and aggregation of particles under the 
combine action of a constant and a gradient external magnetic held acting si- 




Computational Modeling of an MRI Guided Drug Delivery System. 


9 



Simulation Experiment Simulation Experiment 


(a) 


(b) 



(c) (d) 


Fig. 6: Box plots show variations in the major axis length of aggregates from the 
present simulations and the measurements of Ref. [^. The boxes captures the 
lower quartile, median (red line) and upper quartile values. The black line that 
extended from the box show the upper and the lower fence. The blue dots shows the 
average length of the aggregations in each case, a) 0.563 mg/ml^ b) 1.125 mg/ml^ 
c) 2.25 mg/ml^ and d) 4.5 mg/ml. 


multaneously. The present results were compared against the experimental mea¬ 
surements and simulations of for the motion of a magnetic suspension of 
polystyrene magnetic particles (5.5 /xm) with density of 1050 kg/vr? and concen¬ 
tration of 25 kgjw? in distilled water. 

Driven by the uniform magnetic field and the magnetic gradient along the 
a:—axis, the microparticles are successively aggregated and simultaneously trans¬ 
ported in the direction of the magnetic gradient. Four instances of the particle 
locations that show their motion and the mechanism of the formation of the ag¬ 
gregates are presented in Fig. Depending of their close neighbors positions, 
particles may approach others or stay isolated. The duration for the formation of 
the aggregates is t = 1.25s as depicted in Fig. and for larger times only transla¬ 
tional motion of the aggregates is found. Thus, as it is expected, it is found that 
longer time is needed for the formation of aggregates than in the first case studied 
due to the slower magnetic forces from the weaker horizontal constant magnetic 
field. The results of the present simulation are summarized and compared against 
the experimental data and simulations of Ref. in Table 

It is observed that the present results show good qualitative and quantitative 
agreement mainly in terms of the mean velocity and the mean length of aggrega- 
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Table 2: Summary of the results from Case 1 



Particles per simulation 

311 

385 

336 

307 

Particles per aggrega¬ 
tion 

0.563 
mg! ml 

1.125 
mg ! ml 

2.25 
mg ! ml 

4.5 

mg ! ml 

1 

26 

9 

0 

0 

2 

24 

13 

3 

0 

3 

18 

5 

3 

1 

4 

5 

11 

3 

1 

5 

7 

5 

3 

0 

6 

7 

6 

2 

0 

7 

7 

5 

3 

0 

8 

2 

4 

0 

0 

9 

2 

2 

2 

1 

10 

1 

4 

2 

1 

11 

0 

2 

11 

1 

12 

0 

2 

0 

2 

13 

0 

1 

0 

0 

14 

0 

0 

1 

0 

15 

0 

0 

0 

1 

16 

0 

0 

2 

0 

17 

0 

1 

0 

1 

18 

0 

0 

2 

0 

19 

0 

1 

0 

0 

21 

0 

0 

1 

0 

23 

0 

0 

1 

1 

25 

0 

0 

1 

0 

26 

0 

0 

0 

2 

27 

0 

0 

0 

1 

31 

0 

0 

0 

1 

32 

0 

0 

1 

1 

35 

0 

0 

0 

1 


Table 3: Comparison of the present results against the experimental measurements 
from Ref. [^. 


Concentr. {mglml) 

Particle number 

Mean length {gm) 

std 

Mean length {gm) 

std 

- 

- 

simulation 

sim. 

experiment 

exp. 

0.563 

311 

34.02 

23.72 

31 

16 

1.125 

385 

59.14 

43.19 

59 

36 

2.25 

336 

108.93 

83.79 

137 

85 

4.5 

307 

201.49 

109.42 

317 

195 


tions in comparison to the experimental one. The small number of particles that 
is used in the experiment and the simulation is probably the source of the small 
difference in the comparison because initial locations of particles may be crucial 
for the formation and motion of chains. Moreover, results may be inffuenced also 
by the diameter of the particles that was kept constant in the present simula¬ 
tion, while in the experiment had a normal distribution. Consequently, existing 
small particles were easier and faster concentrated around bigger particles and 
thus bigger chains of aggregations may be formulated. 
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Table 4: Comparison of experimental and simulation results for Case 2 


Case 

Mean size 
(particles) 

std size 
(particles) 

Mean velocity 
(pm/s) 

std velocity 
(pm/s) 

Experiment, Ref. |13| 

7 

5 

7.5 

1 

Present 

5.26 

1.89 

8.3 

1.4 

Numerical, Ref. Il3| 

10 

4 

9 

2 



0.0002 


0.00015 


(/) 

X 

(0 


0.0001 


5e-05 


-1-1-1-1-1-1- 

p: ^ Initial positions 

- ^ time = 0.51s 

1-1 

o 

o 

1- 

\ <9 o 

% 



j : 6^ 
csP ■. 

6> 

c 

) 

] . o ;o°Q ;o 

0 oo°° :°°° 

a °0 :0 : 

■ ■ D i 

°Si 

k 

o 

j 

) 


.8 . 

0 

o 

. °i 

1 


i o i 

_1_1_1_1_1_1 

1_1 

i _ i 

i_ 


2e-05 


0.0001 

X axis 


0.0002 


(b) 



Fig. 7: Staring positions of particles (red circles) and final positions (green circles) 
for particles displacement between t = 0 s and: a) t = 0.12 s, b) t = 0.51 s, c) 
t = 0.94 s, and d) t = 1.25 s. 

















































































































12 


N.K. Lampropoulos, E.G. Karvelas, I.E. Sarris 


4 DISCUSSION 

In this study, we were able to validate a new numerical model for the MRI guided 
drug delivery system that is based on magnetic nanoparticle aggregations for the 
navigation of paramagnetic nanocupsules. The present method is based on com¬ 
putational fluid dynamics and it is as general as possible and thus can be easily 
extended for use under realistic in-vitro or in-vivo conditions. The purpose of this 
work is to examine under which conditions the model is reliably working. For 
this reason, two basic particle suspension flows were selected where particles are 
influenced by a constant and in the second case a gradient magnetic held. 

Simulations with the present method is more expensive as the number of par¬ 
ticles increases because magnetic moments should be calculated for every particle 
and for every time step. However, given a particle suspension and the magnitude 
of the external magnetic held the method can be applied and provide very good 
accuracy for the prediction of the size of conglomerates. As it is observed from 
Table when the ratio of the number of the particles to the suspension con¬ 
centration is high, numerical results are very close to the experimental. As this 
ratio decreases, the method underestimates the size of the aggregations but still 
is capable to predict qualitatively their growth. 

An important advantage of this work is also that the developed model is found 
to predict more accurate the dynamics of the aggregates than other simulations 
with similar models as Table presents. Results show that the most important 
quantity for the driving of the nanoparticles, i.e. the mean velocity which is devel¬ 
oped under gradient magnetic flelds, is very closely predicted to the experimental 
values. This improved feature of the proposed model is due to the evaluation of the 
more realistic aerodynamic coefficient C^A^ff that incorporates the effect in the 
drug from the exposed shape of the aggregates in a better way, since it is general 
true that C^A^ff < C^A. 


5 CONCLUSION 

A numerical model for the simulation of magnetically guided drug delivery was 
developed in order to predict the motion of magnetic particles for medical ap¬ 
plications. This method aims at predicting the aggregation’s procedure and the 
velocity of each particle inside the arteries and arterioles of a human body. It aims 
also at simulating the motion of aggregations when these are driven by MRI mag¬ 
netic gradient coils. A novel method is developed for the calculation of the drag 
coefficient that takes into account the exposed area of the particle in the fluid. 

The method was tested through comparison against experimental and numer¬ 
ical data. It was found that the present method can simulate satisfactory the 
experiments in a stationary fluid under a steady magnetic field. Furthermore, the 
model was tested for the acceleration of aggregated particles under the influence 
of a constant and a superimposed gradient magnetic field. The results were very 
close to existing experimental data in terms of velocity and aggregation size and 
comparable to the results from existing simulations. 
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